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ABSTRACT 


The  specific  impetus  for  this  work  was  a  conceptual 
design  of  a  submarine  using  the  toroid  as  the  pressure  hull. 
The  designers  could  not  find  a  ready  body  of  knowledge  to 
obtain  scantlings  for  thier  pressure  hull. 

This  work  began  with  a  review  of  efforts  to  solve 
complete  toroidal  structures.  Several  works  were  found  which 
addressed  general  shells  and  extended  into  partial  toroids, 
but  the  solution  of  a  complete  toroid  was  not  found  to  be  a 
common  exercise.  Some  of  the  these  works  are  briefly 
reviewed.  An  attempt  was  then  made  to  solve  for  the 
displacements  in  a  thin  walled  circular  toroid  using  the 
energy  method.  Several  problems  were  identified  associated 
with  the  structure  geometry  which  make  the  solution  for  the 
complete  toroid  difficult.  In  addition,  the  functional  used 
for  the  energy  method  needs  to  be  more  complex  than  the 
simple  trigonometric  or  power  series  functionals  used  in  this 
work.  These  two  areas,  geometry  and  functionals,  are  fertile 
areas  for  further  study  _ 
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CHAPTER  i 


INTRODUCTION 

The  purpose  of  this  work  was  to  determine  the  response 
of  a  thin  toroidal  shell  structure.  The  initial  goal  seemed 
achievable  until  the  literature  was  reviewed.  A  simple  and 
understandable  solution  was  not  available.  The  disparity  in 
works  was  apparent  in  such  basic  elements  as  coordinate 
systems.  The  solutions  were  generally  directed  towards 
displacements  or  forces.  The  solution  paths  varied  from 
complex  algebra  to  finite  difference  method.  The  solutions 
were  as  varried  as  the  number  of  individuals  attempting  them. 
It.  was  impossible  to  expect  any  simple  correlation  among  the 
few  people  who  attempted  to  address  this  topic. 

The  result  of  this  work  is  the  identification  of  two 
areas  where  more  time  and  effort  is  required. 

A.  The  impact  of  geometric  curvature  on  the  solution. 

B.  The  impact  of  the  shell  geometry  on  the  assumed 
functional  used  in  the  energy  method  of  the  solution. 

To  be  very  specific,  these  two  areas  have  precluded  this 


effort  from  achieving  the  initial  purpose  of  this-  year  long 


effort. 

A  word  needs  to  be  addressed  towards  assumptions  that 
are  made.  The  assumptions  are  fairly  standard  for  a  shell 
being  analysed  in  the  linear  elastic  range.  The  assumptions 
are: 

1.  The  shell  is  made  of  isotropic  and  homogenious 
material  which  obeys  Hook's  Law. 

2.  The  thickness  of  the  shell  is  constant. 

3  The  thickness  of  the  shell  is  small  compared  to 
the  radii  of  curvature. 

4.  A  straight  line  normal  to  the  middle  surface 
before  deformation  remains  straight  and  normal  to  the 
middle  surface  after  deformation  and  retains  its 
original  length. 

In  order  to  establish  a  consistent  presentation, 
geometry  will  be  discussed  here  and,  unless  otherwise  noted, 
all  representations  of  geometry  will  follow  this  coordinate 
system  and  the  definitions.  Solutions  by  other  individuals 
will  be  translated  into  this  reference  system. 
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DEFINITIONS 


a  Ratio  of  R  to  r:  a  = 

Since  R  >  r  to  form  a  complete  toroid,  a  will  always 
be  greater  than  1 . 

C  Circumference  of  the  parallel  circle:  C  -  2  n  r 

r  Radius  of  curvature  in  the  0  direction  of  the 

unloaded  shell  envelope,  or  radius  of  the  circle 
which  is  rotated  about  the  Z  axis  Cat  distance  R)  to 
form  the  toroid. 

tic 

r  :  Radius  to  any  point  on  the  unloaded  toroid  measured 

perpendicularly  to  the  Z  axis. 

r*  =  R  +  r  sin©  =  r  (a  +  sin©). 

R  Radius  of'  rotation  about  the  Z  axis  of  the  circle  of 

radius  r  used  to  form  the  toroid. 

t  :  Thickness  of  the  plate 

c>  Angle  of  rotation  about  the  Z  axis  measured  counter 

clockwise  from  the  positive  X  axis  in  the  X-Y  plane. 

©  Angle  of  rotation  measured  clockwise  from  the  local 

perpendicular  to  the  X-Y  plane  in  the  positive  Z 
director  at  distance  R  from  the  origin  in  any 
direction  <?.  See  Figure  1.  It  sounds  like  a 


lot.  but  it  is  not. 


1  !»'  &*'  I.'  «-■  1.1  l.«  L<  i.ij 
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P  :  Radius  of  curvature  in  the  <p  direction. 

p 

P  =  — —  +  r 
sin© 

Coordinate  System:  The  standard  right-handed  XYZ  coordinate 
system  is  the  global  reference  system.  Other 
coordinates  (eg.  0  and  0)  build  on  this  basic 
reference  system. 

Geometric  Curvature:  Curvature  of  the  shell  due  solely  to  the 
geometry  or  curvature  of  the  unloaded  toroid. 

Load  Curvature:  Curvature  of  the  shell  caused  by 

displacements  in  response  to  the  applied  load.  This 
could  also  be  looked  upon  as  the  change  in  curvature 
from  the  reference  (unloaded)  curved  surface. 

Axis  of  Symmetry:  Axis  about  which  the  small  circle  is 
rotated  to  form  the  toroid  -  the  Z  axis. 

Plane  of  Symmetry:  Plane  which  bisects  the  toroid  -  the  X-Y 
plane 

Meridian:  Intersection  of  any  plane  containing  the  Z  axis 

with  the  toroid.  On  either  side  of  the  Z  axis,  the 
intersection  results  in  a  circle  (for  the  unloaded 
toroid)  with  ©  as  the  variable.  (See  Figure  2) 
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CHAPTER  2 


HISTORICAL  VIEW  AND  THE  PROBLEM 

Historical  View: 

Upon  starting  to  research  the  analysis  of  circular 
toroidal  shells,  the  field  was  found  to  be  not  very  large, 
books  on  shell  theory  and  analysis  would  neglect  or  only 
briefly  mention  the  complete  toroid.  This  could  imply  that 
the  toroid  was: 

< 1 )  a  simple  extension  of  general  shell  theory  or 

C2)  much  more  complex  than  general  shell  theory. 

The  latter  seems  to  be  the  case. 

Senjanovic  11]  gives  a  brief  history  of  shell  theory 
in  his  book.  His  Russian  pride  shows  in  some  instances 
C ". . . shortcomings  were  eleminated  by  representatives  of  the 
Russian  School..."),  but  the  relative  youthfulness  of  the 
field  is  obvious.  Initial  work  in  what  is  now  called  shell 
theory  was  by  H.  Aron  Cl 874)  and  A.  Love  Cl 888)  a  mere 
century  ago.  Senjanovic  does  cite  Flugge  -  another  work 
referenced  here  -  as  one  of  the  developers  of  the  classical 


Cor  western)  theory  of  shells.  Yet  even  Flugge  devotes  only  four 
pages  12]  to  the  toroid.  The  following  works  are  illustrative  of 
the  background  on  toroidal  shells. 

V.  Flugge:  Flugge  starts  with  a  generalized  shell  and 
establishes  force  equlibrium.  From  this  equlibrium 
he  solves  differential  equations  to  obtain  solutions 
for  forces  and  displacements.  His  discussion  of  the 
toroid  points  out  some  of  the  problems  inherent  in 
the  full  toroidal  solution. 

L.  Sobel:  Sobel  131  was  an  understudy  of  Flugge. 

His  doctorial  dissertation  starts  with  the  same 
generalized  shell  theory  but  is  slanted  towards  the 
buckling  analysis  of  the  toroid.  His  solution  was 
worked  on  a  computer  and  is  difficult  to  follow. 

I.  Senjanovic:  Senjanovic  starts  with  a  force 
ballance  on  a  differential  element  in  a  localy 
orthogonal  coordinate  system,  then  uses  several 
methods  to  solve  them.  Most  notably,  he  employed 
complex  variables  and  some  computer  solutions  circa 
1970. 

S'.  Timoshenko:  Timoshenko  141  starts  with  a  plate 


and  develops  into  a  generalized  shell.  His  attention 
to  the  toroid  is  of  the  same  order  as  Flugge.  This 
writer  found  Timoshenko's  general  development  of 
theory  easier  to  follow  than  others  and  hence  the 
references  to  Timoshenko  in  this  area  may  be  more 
frequent  than  to  other  authors  listed  here. 

E.  Tsui:  Tsui  t5]  starts  with  the  constitutive 
relationships  and  develops  them  into  forces.  However 
for  the  toridal  shape  his  solutions  are  on  globally 
orientated  forces  and  displacements.  Tsui  intended 
his  work  more  as  a  handbook  than  a  textbook.  His 
solutions  are  tables  of  values  from  computer  analysis 
of  toroidal  shells.  He  does  provide  however,  a 
consistant  development,  though  shorter  and  in  less 
detail  than  some  others,  for  various  shells  of 


revolution. 


The  Problem: 


The  specific  problem  to  be  addressed  in  this  thesis  is 
the  response  of  a  toroidal  shell  to  hydrostatic  presure. 

More  specifics  on  the  loading  will  be  addressed  in  chapter 
three.  This  is  to  set  the  stage  for  the  problem  itself 
and  lay  out  the  difficult  areas. 

The  goal  of  any  structural  solution  is  to  take  into 
account  the  loading,  geometry  of  the  structure,  and  the 
material  in  the  structure  and  determine  the  response  of  the 
entire  structure.  The  response  is  then  compared  to  a  desired 
standard  or  a  set  of  failure  criteria  and  adequacy  of  the 
design  is  determined.  This  approach  is  generally  slanted 
towards  displacements  of  the  structure  since  only  through 
displacements  are  the  applied  loads  distributed  and 
equi 1 i bra ted. 

In  any  structure,  the  various  elements  are  in 
communication  with  one  another.  The  beam  under  a  simple 
tensile  axial  load  is  the  easiest  example,  but  a  plate  or 
shell  must  also  satisfy  this  connectivity.  Unfortunately, 
the  complexity  of  this  connectivity  or  communication  among 
elements  increases  rapidly.  In  the  axially  loaded  (tensile) 


beam,  only  the  axial  displacement  Csay  the  X  direction)  is 
important.  Load  the  same  beam  laterally  and  now  the  X  and  Z 
directions  are  required  to  describe  the  beams  response. 

Going  to  a  plate  now  requires  the  X,  Y,  and  Z  directions. 

The  shell  now  takes  the  plate  out  of  the  convenient  XY  plane 
for  a  reference. 

With  appropriate  coordinate  system  selection,  some  very 
practical  shells  can  be  handled  easily.  The  cylinder  and  the 
sphere  are  classical  examples.  The  transition  to  a  toroid 
would  -  by  simple  implication  -  be  very  easy.  Start  with  a 
cylinder'  of  radius  r  and  length  of  2nR  and  bend  it  around 
upon  itself.  Discounting  a  few  wrinkles  on  the  inner 
sections  of  the  cylinder,  we  now  have  a  toroid.  The  wrinkles 
however  carry  a  much  deeper  implication. 

The  structure’s  ability  to  distribute  the  applied  load 
changes  drastically.  The  complete  toroid  has  no  boundary 
conditions  comparable  to  the  cylinder's  boundary  conditions 
at  X*0  and  X»L.  These  are  replaced  by  connectivity  Cor 
continuity)  requirements.  The  analysis  of  toroidal  segments 
is  facilited  by  the  ability  to  insert  boundary  conditions. 

But  even  this  facility  can  be  of  limited  use  if  too  much  of 


the  structure  is  included  in  the  analysis.  The  simple 
thin  walled  analysis  on  a  cylinder  for  a  force  balance 
becomes  an  algebraic  exercise  of  no  small  feat  for  the 
toroid.  (This  will  be  shown  in  chapter  five.) 

It  should  be  noted  that  all  of  the  individuals  above 
stated  or  implied  some  difficulty  in  obtaining  a  solution  for 
a  complete  toroid.  The  most  blatant  was  Tsui  E6J  when  he 
states  that  in  problems  involving  the  complete  toroid,  the 
crowns  (  ©  ■  0°and  180°  -  see  FIG  1)  can  not  be  solved 

due  to  singularitj es.  Flugge  171  gives  a  solution  for  the 
forces  but  states  that  they  cannot  be  used  "...  in  the 
vicinity  of  the  top  and  bottom  circles"  (Tsui's  crowns) 
"without  additional  bending. . . ".  He  then  provides  a 
displacement  expression  containing  terms  which  blow  up  at 
©  -  0°  and  180°. 

With  all  of  this  going  against  an  easy  solution,  two 
courses  of  action  were  undertaken: 

A.  The  singularities  at  the  crowns  could  be  avoided 

by  using  some  expression  which  avoids  the  offensive 

term  ( — - — )  (see  Chapter  6).  This  was  avoided  by 
sin© 

referencing  the  differential  element  to: 


-V*  fMinim*  -  J  -  *  - J  *jnui  TJI  njiTJ  «J  ■.«  'jl'J  'J  ’JT-J  WWW  VTT  wy,  rr.1^  IP.'^TWl'r 


B 


r  *  fc  +  r  sin  ©  ■  r  C  a  ♦  sin  ©  > 

The  biblical  addage  still  holds:  If  thine  eye  offends 
thee,  pluck  it  out. 

B.  The  method  of  solution  selected  would  be  one  to 
avoid  further  problems  as  pointed  out  by  Flugge  t8J. 
The  energy  method  would  avoid  messy  mathematical 
stumbling  blocks  in  the  solution,  since  the  form  of 
the  solution  is  assumed  going  into  the  energy  method. 
This  was  almost  true  as  shall  be  pointed  out  latter. 
Additionally,  this  writer  has  not  seen  this  attempted 
by  others.  Therefore,  when  the  line  forms  to  say 
that  this  method  will  not  work,  it  should  be  a  very 
short  line  indeed. 
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CHAPTER  3 


CONSTITUTIVE  RELATIONSHIPS 

An  expression,  or  series  of  expressions,  relating  the 
material  motion  to  properties  such  as  stress  and  strain  Cand 
further  to  forces  and  moments)  is  required.  This  follows  the 
definition  of  constitutive  equation  of  the  material  as  given 
by  Fung  19].  "...mathematical  expression  of  the  mechanical 
property  of  a  material." 

First  we  must  define  our  displacements  (see  Figure  3). 
Three  displacements  will  be  used: 

u  : Displacement  in  the  positive  ©  direction 
v  :  Displacement  in  the  positive  <p  direction 
w  : Displacement  in  the  direction  of  the  positive 
(ie.  outward  pointing)  normal  at  the  point  of 
interest. 

Rotations  will  be  defined  as  follows  (see  Figure  3): 

:  Rotation  of  a  unit  norma]  about  the  tangent  to 
the  meridian  in  the  positive  u  direction:  that  is. 
rotation  about  the  positive  ©  axis. 


cj  :  Rotation  of  a  unit  normal  about  the  tangent  to 

V 

the  parallel  circle  in  the  positive  v  direction; 
that  is.  about  the  local  positive  0  axis, 
o.  :  Rotation  about  the  positive  w  direction. 

The  loading  will  impact  greatly  on  the  final  outcome. 

For  this  problem,  the  loading  will  be  assumed  to  be  a 
constant  hydrostatic  pressure.  This  has  several 
implications; 

A.  Hydrostatic  loading  always  presents  a  loading 
normal  to  the  surface. 

B.  Because  of  the  global  nature  of  the  hydrostatic 
loading,  any  benifit  to  be  garnered  from  symmetric 
conditions  would  be  applicable  to  this  problem. 

The  second  aspect  of  the  loading,  axisymmetric  loading, 
results  in  the  number  of  unknowns  being  reduced  dramatically. 
Both  Timoshenko  1101  and  Flugge  Ill!  state  that  the  total 
deformation  can  be  represented  by  only  two  displacements;  one 
normal  to  the  surface  and  one  along  the  meridian  in  the 
positive  0  direction.  Flugge  also  states  [121  that  that 
stresses  are  independant  of  <2>  and  in  fact  anv  derivatives 
with  respect  to  0  are  zero.  Timoshenko  1131  agrees,  but 
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the  derivative  going  to  zero  is  implied,  not  stated  as 
bluntly  as  in  Flugge. 


To  take  this  one  step  further,  a  lot  of  things  will  go 

to  zero  because  of  axisymmetric  loading.  Sheer  stress  and 

strain  are  cross  coordinate  terms,  ie. 

[*]A  =£_(£_[*]  ) 

6<t>  66  6<p 

but  aill  ^  50  l*300  =  l*]*e  =  0. 

Because  of  the  type  of  loading  described,  v  and  q,  are  not 
required  for  the  description  of  motion.  v  is  not  required 
because  of  the  independance  of  the  motion  on  4>,  and  o.  is  not 
required  because  it  depends  on  a  partial  of  <p.  All  of  this 


is  the  result  of  axisymmetric  loading. 


STRAINS 

Midplane  strains  are  as  follows: 

Cq  :  Is  the  strain  in  the  positive  0  direction.  Cq 
is  made  up  of  two  parts  (see  Figure  4): 

1 .  Cq  due  to  u 

At  two  points  on  the  meridian:  P  and  Q 

Up  =  U  uQ  =  u  +  is  de 

€  =  —  =  — —  Cu  +  —  de  -  u)  =  -  — 

1  r  d0  66  r  66 

2.  Cq  due  to  w 
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r*vw 


r#. t 


.  Ala  Atm  *1-  »««  |K  |V>|  < 


‘il.’iL  it 


Length  of  a  segment  is  equal  to  the- radius 


times  the  angular  rotation.  The  original 


length  is  r  d0.  At  point  P  the  radius 


increases  by  w  to  Cr  +  *0  so  the  new  length 
would  be  Cr  +  w)  d0.  Similarly  at  point  Q  the 


radius  increases  but  is  now 


Cr  +  w  +  —  d9).  The  new  length  would  be 
<50 

Cr  +  w  +  £J2d0>d0.  Averaging  the  two  lengths 
66 

at  P  and  Q  gives  an  average  A1  =  wd0  +  ^ 
Dividing  by  the  original  1  =  r  d0  gives: 


c  =  h  Iw  +  i.  ~  d01  Neglecting  second 
r  2  <50 

order  terms  yields:  c  =  " 


1  <5w 


Total  Cq  is  therefore: 


=  rr  t  — —  +  W 1 


'0  "  F 


Cp  :  Is  the  strain  in  the  positive  <p  direction 
Csee  Figure  4).  Since  the  body  is 


symmetrically  loaded,  all  0’s  have  the  same 
displacements.  Therefore  can  be  calculated 
by  calculating  the  change  in  the  circumference 
of  the  parallel  circle  at  <p  -  constant.  Since 


the  circumference  is  related  to  the  radius. 
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the  radius  (r  >  will  be  used  to  determine 


strain  c^. 


C  =  2  n  r 


* 


At  point  F*.  displacements  w  and  u  contribute 


>H 

to  the  Ar  as  follows: 

Ar  =  u  cos©  +  w  sin© 
AC  =  2  n  Ar* 


c  ^  =  —  =  i-jp  Lu  cos©  +  w  sin©] 
C  r 


ROTATIONS’ 

is  determined  from  two  points,  P  and  Q,  on  a  meridian 

V 

seperated  by  displacement  u  in  the  positive  ©  direction  Csee 
Figure  5).  The  unit  norma]  vector  at  P  must  rotate 
through  d©  to  go  from  P  to  0  due  to  displacement  u,  ie. 
r  d©  =  u 


If  w  was  constant  from  P  to  Q,  no  additional  rotations  would 
be  required  since  the  surface  displaces  uniformly  through  d©. 
If  *  0.  then  a  rotation  of  the  local  normal  at  0  will  be 
required  to  maintain  its  normal  relationship.  This  rotation 
w  i  ] ]  be • 


£j2d©/(  r  d©)  = 


1  6w 


This  makes  the  total  • 

<P 


_  1  f 6\t>. 

0  r  66 

Because  of  the  symmetric  loading  and  no  subsequent 

variations  in  the  displacements  with  <p,  no  v  displacements  or 
<5w 

~  exists.  However  a  can  exist  (see  Figure  6). 

66  ^ 

Due  to  symmetry  in  loading,  Iw^J  is  constant  in  <p. 
However,  for  to  change  direction,  an  additional  rotation 
perpendicular  to  must  exist.  This  is  the  rotation 

ci0=  -  cos0  d <t> 

Figure  6  gives  a  graphical  presentation  of  Oq. 

Timoshenko  L141  gives  a  similar  result  but  achieves  this  as 
a  side  result  of  obtaining  curvature. 


LOAD  CURVATUtLS 

Load  curvatures  are  as  follows: 

Mq.  Change  in  curvature  along  the  positive  0  axis 
due  to  the  applied  load. 

x^.  Change  in  curvature  along  the  positive  <p  axis 


due  to  the  applied  load. 


Figure  6 


Rotation  Qq 


A  u)  -  Sin  (J $ )  >=  u)t 


It  is  significant,  to  this  writer  at  least,  to  seperate 
these  curvatures  from  the  geometric  curvature  defined  in 
Chapter  1.  The  geometric  curvature  gives  one  an  indication 
of  how  difficult  the  problem  will  be  L 15 1  and  defines  the 
initial  or  unloaded  reference  surface  for  further  work.  The 
load  curvature  is  the  sole  result  of  the  structure’s  response 
to  the  load  applied  and  directly  related  to  the  displacements 
use  to  equilibrate  the  load. 


parameters  in  terms  of  displacements.  Adding  to  this  list 
the  general  expression  for  planer  stress,  we  get  the 


following: 


i£< 
/ 1 


8 


fl 


k'JM 


I 


a 


& 


$ 

-J 


_  (ee  + 


ci  [£g  ♦  w]  +  - - - 

1-r*2  r  rCa  +  sin©) 


[u  cos©  +  w  sin©]) 


o,  =  -L 


[u  cos©  +  w  sin©]  +  —  [—  +  *f)> 


1-v  r(a  +  sin©) 


r  <5© 


This  almost  completes  our  constitutive  relationships. 


Stress,  strain,  rotation,  and  local  curvature  can  be 


expressed  in  terms  of  the  displacements  u  and  w.  The  last 


step  is  to  proceed  to  normal  forces  and  moments.  This  is  as 


f ol lows: 


o  dz 

t 


-  e,  * 


o  z  dz 
\ 
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CHAPTER  4 

EQULIBRIUM  OF  FORCES 

The  next  aspect  of  the  solution  to  be  presented  will  be 
that  of  establishing  equlibrium  conditions  for  the 
differential  element.  The  differential  element  is  shown  in 
Figure  7.  Geometry  rears  it's  ugly  head  in  that  the 
element,  needs  two  reference  points  to  describe  it’s  surface. 
The  first  is  the  axis  of  rotation  of  the  shell  in  the  B 
direction.  The  second  is  the  origin  of  the  global  coordinate 
system.  The  first  axis  can  be  referenced  to  the  global 
coordinate  system  by  R  and  4>. 

Solutions  for  partial  toroidal  shapes  all  take  advantage 
of  the  known  boundary  conditions  on  the  edges.  When  a 
complete  toroid  is  analyzed,  the  boundary  conditions  double 
back  upon  themselves  and  become  internal  vice  external 
forces.  In  these  instances,  symmetric  loading  and 
deformations  are  invoked.  Tsui  L17J  gives  tabulated 
computer  results  for  both  solutions  and  the  differences  are 
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FI**** 


Differential  Element 


obvious . 


In  the  following  discussion,  and  the  associated  figures, 
the  following  conventions  will  be  used  Csee  Figure  8): 

9 

F  is  a  generalized  force  ie.  a  force  or  moment  in 
the  i  direction  <i  ■  ©,  tf>,  n>.  The  units  for  Fv  are 
t Force/unit  length!. 

F  is  the  force,  ie. 

I  7 

F  *  [Force/unit  length]*! length!  ■  [Force! 

Ne-  *V  Normal  forces  in  the  <^>  direction  acting  on 
the  side  of  the  element  on  which  <^>  is  constant. 

Qe  Sheer  force  acting  on  the  4>  ■  constant  face  of 
the  element.  Note:  Due  to  symmetry  of  loading, 

%  * 

0 

M0'  Moment  to  induce  rotation  in  the  <^> 

direction.  The  moment  is  defined  as  positive  with 
respect  to  the  right  hand  rule  as  applied  about  the 
local  positive  <^>  direction. 

Each  force  will  be  addressed  graphically  with  three 
Isometric  diagrams  and  a  geometric  breakdown  to  the 
appropriate  F  .  The  F  for  each  force  CNe,  N  .  . ect. )  can 

then  be  summed  for  a  force  balance. 

-  33  - 


Figure  10 


Contribution  to  F  doe  to  N>: 

v  <p 


F^:  -N^  +  cosa  cos/3 


F0:  sina  cos/3  -  cos/3  sina 


F  *.  N  .  sin/3  -  N  ,  sin/3 


Figure  11 


Contribution  to  Fv  doe  to  Q^: 

9  9  9  9 

F^:  coso  sin/3  ♦  coso  sin/3 

9  9  9  9 

F©:  Q©  sina  cos/3  ”  Q©  cosa  sin/3 

>  9  9  9 


F.J  Qa  cosa  cos/3  ♦  cosa  cos/3 


Figure  12 


P 


Contribution  to  F^  doe  to  Z: 

F>:  -  Z  cosa  sin/3 
<P 

F^:  -  Z  cos/3  sina 


-  Z  cos/3  cosa 


Figure  14 


Contribution  to  F  doe  to 

V  ? 

*  >  •  * 
Nwk  cos  ft  cosa  -  CM.  +  — 5C  d©>  cos/3  cosa 

, 

F0:  -  M,  cosft  sina  +  +  — £  de>  cosp  sina 

^  ^  66 

»  >  >  » 

F„:  -  M .  cosa  sin/3  ♦  <M.  ♦  — 2  d©>  cosa  sinf? 


E 


I 


Relationships  between  the  N^and  N  . must  now  be 


established. 


Ne  -  N0  r  dO 


6N, 


N„  +  —12  =  CN'  +  S^Xr*  +  dtp 

60  a  60  *e 


N^r*d0  +  N^^d*  +  (^!2;>r*d0  +  <f^£)c2£. 


’o' 


*6 


60 


60  60 


6N„ 


=  (N@  +  — X)  r  d$>,  neglecting  terms  with 


6© 


being  relatively  small 


Similarly 


N  .  =  N  ,  r  dO 
tp  <P 


'O 


Q  e  r  d  tp 


6©0  *  ^ 
q  +  < — x;  =  <q  +  — X)  rm  dtp 
w  60  w  60 


M0  ■  M0  r  d tp 


M. 


6Ma  .  Mo  k 

— X  e  <M  +  — X)  r  dtp 
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CHAPTER  5 


ENERGY  METHOD 

Any  attempted  solution  must  maintain  some  sort  of 
equlibrium.  In  the  energy  method,  the  equlibrium  comes  from 
minimizing  the  energy  from  the  external  load  and  the  internal 
energy  of  the  shell.  For  the  external  energy,  a  surface 
integral  of  the  hydrostatic  load  and  the  normal  displacements 
gives  this  value.  The  internal  energy  is  determined  by 
integrating  stress  times  strain  throughout  the  shell  volume. 
The  external  energy  is  subtracted  from  the  internal  energy 
and  the  resulting  functional  is  generally  refered  to  as  n, 
the  total  potential  energy  of  the  structure.  With  V 
representing  the  internal  energy  and  V  representing  the 
external  energy,  the  functional  looks  like: 

n  -  v  -  v 

II  the  value  of  II  is  at  a  minimum  value,  the  structure 
will  be  in  equlibrium.  This  is  the  principle  of  total 
potential  energy  1163.  This  work  will  not  address  buckling 
or  post  buckling  behavior,  only  the  linear  elastic  region. 


•*rwrwmr.  mwwv  wn  itoiwt  vt’sr  '.^  yw  g*  im**  gwvi  w  ir«  y  yjj» y  "ji  "Jin  pj  iu». 


For  this  situation,  equlibrium  is  achieved  when: 

sn  _  6V_  _  S¥  m  0 
6q  6q  6q 

where  is  any  generalized  parameter  used  to  determine  V  and 
V. 

Unfortunately,  the  parameters  used  to  determine  V  and  V 
are  u  and  w  -  which  is  what  we  are  attempting  to  solve.  To 
get  around  this  difficulity,  a  solution  will  be  assumed. 

Then  an  energy  balance  will  be  conducted. 

The  assumption  of  a  solution  is  not  without  risk. 

First,  the  solution  must  meet  known  natural  boundary 
conditions.  Then  the  accuracy  of  the  solution  Is  limited  by 
the  closeness  of  the  assumed  solution  to  the  actual  (and 
unknown!  solution. 

In  the  case  of  the  toroid,  there  are  no  "boundary 
conditions"  as  was  pointed  out  in  chapter  four.  In  place  of 
the  boundary  conditions,  we  have  conditions  of  compatability 
to  limit  the  solution.  As  stated  earlier  ([19]  and  [20]), 
svmmetry  of  deformation  can  be  assumed  for  the  hydrostatic 
loading.  The  plane  of  symmetry  is  the  XY  plane. 
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To  determine  the  internal  energy  we  proceed  as  follows: 

V  ■  /  o  c  dV  ;  v  a  j  -  ©,  4> 

V  V  J  l  J  9  J 


1  /_  <N  e  +  M  x  >  dS 

2  S  v  j  v  j  i  j  v  j 


Recall  that,  due  to  the  type  of  loading, 

cross  terms  ^  j>  equal  zero. 


v  -  I  u  CN©  ce  +  N*  e*  +  Ma  *0  +  M*  V  dS 


Eq 

1  : 

Ne 

e©  " 

k 

Ccg  +  v  e 

.  )  £ 
6 

© 

;  K 

,  t  n 

1-u2 

* 

K 

(Cq  +  2  u> 

e© 

e6  + 

e6) 

Simi larly : 

Eq 

2: 

N* 

e+  = 

k 

ie*  +  2  v 

e* 

e©  + 

cp 

Eq 

3: 

M© 

*©  = 

D 

c*©  +  V 

*© 

i  D 

_  E  h3 
12C1-V2) 

s 

K 

crz)  c*e 

+  2 

v  *e 

*0  + 

*P 

Eq 

4  : 

M6 

*6  = 

k 

crz)  c*6 

+  2 

V  *0 

*©  + 

*P 

»c 


Recall  from  chapter  four: 
Eq  5,6:  eA  ■  i  (—  +  w)  ;  e. 


60 


'6 


<u  cos©  +  n  sin©> 


6u  6W 


Eq  7.8:  *e  '  ^  - 


■> ; 


r(a  +  sin©> 

_  cos© 


cu  -  *2) 


r"  6©  60"  *  r2(a  ♦  sin©>  6© 

The  unknowns  u  and  w  still  persist.  The  assumption  will 


be  made  that  u  and  w  are  functions  of  ©.  A  series  solution 
will  be  assumed  as  follows: 
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rwji  Tw^nmrvnf^rm^v.  yrnTTir  w^^’.iWfffiwwrjw.^L^w[ia 


w< ©>  *  I  w  sin(m©>,  limits  of  summation  will  be 

TTi 

addressed  later 

u<©)  ■  £  u  —  <sin(m©>> 
m  60 

■  Z  u  m  cos(m6) 

m 

This  turns  out  to  be  not  a  very  good  assumption  for  the 
the  series  representation  of  the  displacements.  There  are 
some  inherent  limitations;  lor  instance  sincm©)  is  always 
zero  at  ©  *  0.  While  not  in  conflict  with  the  conditions  of 
compatibility,  this  is  an  additional  constraint  on  the 
solution  displacements.  If  the  actual  displacements  are  not 
zero  at  ©  «  0.  then  there  is  a  guaranteed  error  in  the 
assumed  solution.  The  question  to  be  answered,  as  with  all 
assumed  solutions,  is  how  good  is  the  answer  obtained 
compared  to  the  answer  that  is  required.  More  on  this 
subject  in  Chapter  6. 

For  ease  of  notation,  let: 

0  *  sin<m©) 

ff» 

> 

<t>  «  m  cosCm©) 

rr, 


0 . 


-m2 


sin<m©) 
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£ 


-t.  A.  ■/»„  A.  V- 


Therefore: 


Eq  9 :  w<  0 )  ■  Z  w  0 

rr  m 

Eq  10:  u(0)  ■  2  u  4> 

m  m 

Boundary  conditions  for  this  problem  are  actually 

compatabi lity  conditions.  At  0  ■  90°  and  270°  the  u(0> 

displacement  should  be  zero  and  — (w<0>>  should  be  zero. 

SO 

Zw  <p  )  m  z  m  £_<  0  )  >  Z  w  0’ 

66  60  m  m  rr'60m  mm 

Since  w  are  merely  coefficients  now  in  the  series  solution 
thev  are  unaffected  by  the  derivative.  To  look  at  this  in 
more  detail: 

A&j 

—  =  w,  cos<0)  +  w„  2  cos(20>  +  w  3  cos(30>  ♦  .  . 
60  1  2  3 

+  w  n  cos(n0) 

r. 

For  conveince.  let  n  *  10  lor  the  trial  solution.  At  0  ■  9 
and  270°.  cos<m0>  *  0  for  m  ■  1,  3,  5,  .  .  and  cosCmO)  ■  ± 
for  m  «  2.  4.  6.  .  .  .  For  the  series  to  be  equal  to  zero 

set  ■  0  for  m  *  2,  4,  6,  .  .  .  Similar  reasoning  gives 
um  ■  0  for  m  »  2,  4,  6,  .  .  .  (This  explains  why  the 

functional  for  u  was  selected  as  the  derivative  of  the 

rr. 

functional  for  w  . )  As  shall  be  seen  later,  the  method  is 

m 

smart,  enough  to  realize  this  limitation.) 

The  elemental  area  (dS)  must  be  investigated  for  the 


toroid.  dS  ■  dl@  dl^  ;  dl@  ■  r  d6  ;  dl^  ■  <R  +'  r  sin©)  d<£ 
*  r  d©  C.R  +  r  sin©)  dtf> 

As  a  check:  S  -  S™f™v  <R  +  r  sin©)  d©  d* 

■4  n2  R  r.  This  agrees  with  the  CRC 
Standard  Mathematical  Tables  [211  for  the  surface  area  of  a 
toroid,  therefore  dS  is  correct. 

The  ultimate  objective  is  to  get  —  •  0.  In  this  case 

q  «=  u  and  w  .  So  for  each  specific  u  -  that  is  u.  -  one 

v  m  m  r  m  k 

gets  an  expressison  such  that: 

sn  m  6V_  _  SV_  m  0 

6uk  <5uk  6uk 

Likewise  a  series  of  m  expressions: 

<M1  _  *V_  _  6V_  =  Q 

6w,  <5w,  6w, 

k  k  k 

For  a  nominal  m  =  10.  this  yields  20  equations.  Fortunately, 
u  and  w  are  20  unknowns. 

it.  m 

Now  the  internal  energy  can  be  formated.  Equations  9 
and  10  are  substituted  into  equations  5  through  8.  These 
equations  are  then  inserted  into  equations  i  through  4.  This 
is  combined  with  the  expression  for  dS  to  give  V,  the 
internal  energy.  To  facilitate  the  future  derivatives  which 
will  be  required,  the  subscripts  m  and  n  will  be  used  to 


indicate  the  seperate  terms  in  a  product,  ie. : 


rv 


V 


n  r2  K  f2Tt[- —  IZCu  u  0  <f>  +  uw  0  0 

o  L  2  m  r,  rrr  m  m  n  ^  m  n 


+  W  U  <p  <f> 

m  n  m  n 


+  W  W  0  0  ) 

m  n  m  n 

►  - - £ZCu  u  cose  +  u  w  0,0  sine 

r  (a  +  sin6; 

♦  w  u  4>  0  cose  +  w  w  0  0^  sine) 

m  n  m  n  m  n  m  n 

►  - - - -zz<u  u  0  0  cos2e 

r2(a  +  sir.e>2  ™  n  m  m 


>  > 

+  u  w  0  0  cose  sine  +  *»  u  0  0  cose  sine 

m  n  m  r<  m  n  m  n 


+  w  w  <p  4>  sin2e; 

m  r»  m  n 


.2  ••  ••  ••  "  1 
♦  c  n— >IZCU  U  0  0  -  U  W  0,0  -  W  U  0  0 

.  ..  4-  m  n  m  m  m  n  m  n  mn  m  r 


w  W  0  0  > 

m  n  m  n 

.  2  i'  cose 


12r4  Ca  +  sine) 


ii  ;  ft 

:zcu  u  00  -  uw  00 

m  n  ^  m  m  m  n  m  n 


it  ;  »•  7 

“WU  00  +  WW  0  0  ) 

mn  m  n  mn  m  n 

+  C  — — >  - C°gi? - (U  U00  -  UW00 

1 2r4  <a  +  sine;2  m  n  *"  w  m  n  "  n 

y  >  r  y 

-  k  u  00  +ww  00)lCa+  sine;  de 

m  n  m  ^  n  m  n  m  r.  J 

XV 

This  is  all  well  and  good  but  actually  -  is  what  we 

6qv 

really  need.  Again,  for  this  problem  qv  **  ufc  and  *k;  and 

k  *  i,  2,  .  .  .  10.  Applying  term  by  term  differentiation  is 

straightforward  but  tricky  and  tedious.  A  couple  of 


examples: 


2_(wu00;-**00 

tr>  n  m  n  mm 
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But  observe: 


C  u  u  <P  0  )  ■  u  0  0,  +  u  0,  0 

m  n  rr.  r,  m  ^k  nTk~r 


2  "m  K  K 


Since  m  and  n  can  both  serve  as  a  counter  Index,  the  two 


terms  may  be  combined.  The  key  to  this  is  that  both 


functionals  (0  )  have  to  be  of  the  same  order  of  the 


derivative.  Observe: 


-  Cu  u  0  0  )  =  U  0  0,  +  U  0,  0 

rr.  n  m  n  m  T  m  T  k  n  T  k  r  n 


"  “m  C*m  *k  +  *k  <0 


Doing  all  of  the  above  and  collection  like  terms  gives  the 


following  two  expressions: 


The  partial  of  V  with  respect  to  ufc. 


— -  =  rr  k  pu  f*nC 2  <fi" 4>''  +  2  V  9°?? - C0’V  +  h") 

.  ..  L  m  O  rn  k  ^  k  m  m  k 


Ca  +  sin©) 


+  .  2  cos.  0'  0' 
Ca  +  sin©)2  "*  k 


+  <HT)C2  *m*k 


"ul”  2  v  cos© 


r<*k*m  +  KK  cose) 


2  cos  © 


Ca  +  sin©) 


0  0  ))  Ca  +  sin©)  d© 


2  ”rrvk 


+Iw  /2n  C 2  - C0"0  sin© 

m  0  k  m  Ca  +  sin©)  k  m 


+  0  0’  cos©)  ♦  2  cos^sine  > 

Ir  -  2  rk'm 


Ca  +  sin©) 


2  v  cos© 


•  •  »  n  » 


+  <r7>c-2  0,0m  -  -  -  C0  0  «•  0  0) 

1  K  rn  ✓  -  »  ^  j  ^  ’v  K  ^  rv  1  k 


Ca  +  sin©) 


2  cos  © 
Ca  +  sin©) 


—  0k0^))Ca  +  sin©)  d©] 
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The  partial  ol  V  with  respect  to  w 


—  =  n  K  pu  /2nC 2  + 

£..  L  m  O  m  k 


k 

2  v> 


6w 


Ca  +  sin©) 


<0  0,  sin©  +  0  0,  cos©> 

m  K  m  k 


,  2  cos©  sin©  . ’ , 
+  - -  0  0 


Ca  +  sin©) 


2  T  nvk 


2  cos  © 


Ca  +  sin©) 


0  0  ) )  Ca  +  sin©)  d© 


2  ^m^k 


+zw  /*"  <2  -  «  ^.«4ae..  0  * 

rr>  O  k  m  ,  _  ^  k  m 


Ca  +  sin©) 


2  sin  © 


Ca  +  sin©) 


2  ~  k  ~m 


+  Clr-r)C20  0  + 


2  v  cos© 


r?' '“^k^m 
2  cos2© 


Ca  +  sin©) 

2  +  sine)  d©] 


Ca  +  sin©)' 

The  external  energy  CV)  should  be  looked  at  next.  The 


energy  is  the  product  ol  the  load  -  hydrostatic  pressure 
CP)  in  this  case  which  is  always  normal  to  the  surface  - 
times  the  displacement  along  the  line  of  action  of  the  force 
-  wC©)  -  in  this  case. 


V  «  /  P  k  (©)  dS 

m 

«  / 2Tr  /2n  P  I  w  0  rCR  +  r  sin©)  d©  d0 

O  O  m  m 


Substitute  R  *  r  a 


*  2  n  P  r2  I  w  /2n  0  Ca  +  sin©)  d© 


LUVj 


And  the  partials  are: 


2  n  P  r2  /2Tr  0,  (a  +  sin6)  d6 

O  k 


2L  =  o 

6u, 


The  series  expressions  can  now  be  rearranged: 


<513  _  6V  <5W 


<5qt  6qt  <5q^ 


6V  _  <5W 


<5q  <5q 


.  q  «  u 


k  ’  k 


This  can  be  represented  in  matrix  format  as  follows: 


Eq  11 


=  CD!  <u>  +  [E]  <w>  -  0  -  — 


Where  C  t)  3  represents  the  coefficients  of  u  in  - 

<5u. 


(see  page  9). 


L E ]  represents  the  coefficients  of  w,  in  - 

*uk 


<  see  page  9) 


(u>  =  u  =  u4 ,  u,,  u 

m  l  z  d 


1W>T  =  «  Hlf  W2,  W3> 


Similarly  in  matrix  format: 


Eq  12 


QL.  «  [FI  <u>  +  C  G3  <w>  -  <P>  «  — 


Where  t  F  3  represents  the  coefficients  of  u  in 


(see  page  10). 


[G]  represents  the  coefficients  of  w  in 


(see  page  10) 
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<p> 


6V 

<5w, 


SV 

<5w. 


+  + 

<Sw, 


►*SL 

6w . 


k  '"1  ”'2  ''"iO 

Solving  equation  11  for  <u>  gives: 


<u>  =  -  [ D1 


- 1 


[El  <w> 


Use  this  repression  in  equation  12: 

< [ FI  CD!"1  [El  +  [G J >  <w>  «  <P> 

Solve  for  <w>: 

<w>  =  i [ F  3  ID]-1  [  E3  +  t  G3  > ” 1  <P> 

If  the  solution  for  w  and  u  are  correct,  then  these 

mm 

can  be  used  to  reconstruct  w<©>  and  uC©>  and  from  these, 
using  the  constitutive  relationships  in  chapter  four. 


► 


stresses  in  the  6  and  $  directions  can  be  calculated.  Roark 
1221  gives  expected  values  which  can  be  used  for  comparison. 

Some  Justification  for  Roark’s  results  should  be 
provided  since  this  is  how  the  solution  will  be  checked.  He 
provides  no  specific  reference  for  these  answers  in  his  book, 


but  working  backwards  one  can  see  that  a  simple  thin  walled 


approximation  analysis  was  performed,  similar  to  that  of  a 


cylinder.  Roark  provides: 


> 


Pr 

2t 


'e 


Pr 

t 


(See  Figure  15) 
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a  +  xSinG 


a  +  sin@ 


( 


is  a  simple  force  balance  in  which  the  pressure  times  the 


internal  cross  sectional  area  equals  the  stress  times  the 
material  area. 


PA 

P  (n  r2) 
°4> 


O.  A 

<p  m 

(.2  n  r>  t 

Pr 

2t 


A  similar  approach  is  used  on  oQ.  The  pressure  acts  on  the 


area  A,  which  in  this  case  is  (2  r)(2  n  R> .  acts  at  0  and 

e  +  160°.  The  thickness  is  constant  at  t  but  the 


circumference  at  6  and  0  +  180°  are  different. 


CX0)  •  2  «  r*  •  2  n  rCa  +  sin0) 

The  metal  area  at  0  is  therefore  AC0>  *  C(0)  t.  The  force 
balance  is: 


PA*  o<0 CC0  )  t  +oC©_>  C<©^>  t 
11  2  2 

with  0,  «  ©  and  0„  ■  0  +  180°.  Substitute  Roark's  o(0)  into 
1  2 


the  above  equation  gives: 


P  4  n  r  R 


Pr 

t 


a  +  ^sin© 

C - - - -)  2  n  rCa  +  sin©  >  t 

a  +  sin©1 

p  a  +  Jsin© 

- —  c - - - -)  2  n  rCa  +  sin©  >t 

t  a  +  sin©^ 

2 


Cancelling  some  terms  shows  the  equality: 


This  must  have  been  how  Roark  arrived  at  these  relationships. 


This  is  how  the  solution  will  be  checked. 

THE  SOLUTION 

The  method  outlined  above  requires  the  integration  of 
several  quantities  which  are  not  readily  found  in  a  table  of 
integrals.  The  integration  was  done  using  Simpsons  Rule 
[231.  It  was  found,  using  known  integrals,  that  50  steps 
with  in  the  interval  0  to  2 n  gave  good  accuracy  <5  to  6 
significant  digits  accuracy).  A  program  was  written  to 
accomplish  the  integration  and  the  subsequent  matrix 
manipulations  to  solve  for  <u>  and  <w>. 

To  check  the  output  from  the  program,  a  spreadsheet  was 
set  up  to  handle-: 

wc  O)  *=  I  w  <(> 

m  m 

u<.©)*  Z  u  <(> 

m  m 

The  spread  sheet  also  calculated  c Oq,  and  o^.  As  a 

check,  was  averaged  and  used  to  calculate  a  pressure  to 
compare  with  the  input  pressure. 

THE  RESULTS 

The  results  are  disheartening.  In  the  sample  output 
provided,  only  one  of  the  stress  values  used  as  a  comparison 


6wC62 

66 

6u(  6) 
66 


I  w  4> 

m  n 

z  11  0 

'  n 


yyy- 

■B.  -  »  1 


came  close  to  the  expected  value.  o^  should  have  been 
constant  and  have  a  value  of  5000  psl.  It  was  not  constant 
but  it  did  manage  to  equlibrate  most  of  the  applied  pressure 
<99.950.  This  is  pretty  good  for  an  approximation.  The 
is  another  story. 

Oq  is  too  low  and  does  not  follow  the  solution  in  Roark 

P  r  a  +  |sine 
°  t  a  +  sin© 

The  following  are  the  default  values  in  the  program: 

GEOMETRY 

Radius  of  rotation  :  8  inches 

Radius  ol  the  circle  :  2  inches 

Thickness  of  the  shell  :  .02  inches 

MATERIAL 

Poisson’s  ratio  .3  Modulous  of  elasticity 

:  3E+07 

SOLUTION 

Number  of  steps  in  the  Simpsons 

integral  :  50  Number  of  terms  to  be 

used  in 

the  series  approximation  :  10 


Pressure 


:-100  psi 


The  program  output  for  w  and  u  are: 


m 

w 

m  m 

U 

1 

m 

-1 . 213117E-03 

m 

-1.137496E-03 

2 

6. 043412E-09 

1 . 491326E-09 

3 

-1 . 034239E-04 

-1 . 151281E-05 

4 

3. 250964E— 09 

2. 026468E-10 

5 

-3. 933767E-05 

-1 . 579169E— 06 

6 

3. 328626E-09 

9. 227051E-11 

? 

-1  . 16143E-05 

-2. 382749E-07 

8 

3. 831 875E-09 

5 . 9864E-1 1 

9 

-2. 034983E— 06 

-2. 535156E-08 

10 

2.  81034E— 09 

2. 817199E-11 

Notice  that  the  even  numbered  terms,  which  are  to  be  set 
to'  zero  as  indicated  earlier,  are  several  orcferes  of 
magnitude  below  most  of  the  odd  values. 


The  reconstructed  output  from  the  program  is  as 


The  spreadsheet  out  put  to  check  the  solution  is: 


P  « 

-100.00 

psi 

-99.90  < 

Sin  Press 

Cpsi) 

v  » 

0.30 

E 

1-v2 

3. 3E+07 

E  * 

3E+07 

R  * 

8.00 

in  a  ■ 

4.00 

r  = 

2.  00 

in 

t  «* 

0.  02 

in 

o 

w 

u 

6u 

°o 

°+ 

60 

60 

<°> 

<.  in) 

(  in  ) 

<  in) 

,in  > 
(racT 

<psl) 

<psi> 

0 

0. 0E+00 

-1 . 8E-03 

-1. 2E-03 

0.0E+00 

-1461 

-4870 

15 

-4. 4E-04 

-1. 4E-03 

-1. IE-03 

4. 2E-04 

-1710 

-4738 

30 

-7. 2E-04 

-8. IE-04 

-9 . 8E-04 

6. 8E-04 

-1946 

-4607 

45 

-9. 0E-04 

-5. 7E-04 

-7. 8E-04 

8. 4E-04 

-2123 

-4405 

00 

-1.0E-03 

-4. 2E-04 

-5. 4E-04 

9. 6E-04 

-2258 

-4248 

75 

-1. IE-03 

-2. 2E-04 

-2. 8E-04 

1 . 0E-03 

-2345 

-4159 

90 

-1 . IE-03 

-2. 3E-08 

-2 . 8E— 08 

1 . IE-03 

-2375 

-4131 

105 

-1 . IE-03 

2. 2E-04 

2. 8E-04 

1 . 0E-03 

-2345 

-4159 

120 

-1 . 0E-03 

4. 2E-04 

5. 4E-04 

9. 6E-04 

-2258 

-4248 

135 

-9. 0E-04 

5. 7E-04 

7. 8E-04 

8. 4E-04 

-2123 

-4405 

150 

-7. 2E-04 

8. IE-04 

9. 8E-04 

6. 8E-04 

-1946 

-4607 

165 

-4. 4E-04 

1 . 4E-03 

1 . IE-03 

4 . 2E-04 

-1710 

-4738 

180 

-9. 6E-08 

1 . 8E-03 

1 . 2E-03 

9. 2E-08 

-1461 

-4870 

195 

4 . 4E-04 

1.4E-03 

1 . IE-03 

-4. 2E-04 

-1269 

-5190 

210 

7. 2E-04 

8. IE-04 

9. 8E-04 

-6. 8E-04 

-1085 

-5498 

225 

9. 0E-04 

5. 7E-04 

7. 8E-04 

-8. 4E-04 

-894 

-5654 

240 

1 . 0E-03 

4. 2E-04 

5 . 4E-04 

-9. 6E-04 

-747 

-5769 

255 

1 . IE-03 

2. 2E-04 

2. 8E-04 

-1 . 0E-03 

-659 

-5854 

270 

1 . IE-03 

6. 8E-08 

8. 4E-08 

-1 . IE-03 

-630 

-5886 

285 

1 . IE-03 

-2. 2E-04 

-2. 8E-04 

-1 . 0E-03 

-659 

-5854 

300 

1 . 0E-03 

-4. 2E-04 

-5. 4E-04 

-9. 6E-04 

-747 

-5769 

315 

9. 0E-04 

-5. 7E-04 

-7. 8E-04 

-8. 4E-04 

-894 

-5654 

330 

7. 2E-04 

-8. IE-04 

-9. 8E-04 

-6.8E-04 

-1085 

-5498 

345 

4. 4E-04 

-1  .  4E-03 

-1 . IE-03 

-4. 2E-04 

-1269 

-5190 

360 

1 . 9E-07 

-1 . 8E-03 

-1 . 2E-03 

-1 . 8E-07 

-1461 

-4870 

Note:  Due  to  the  precision  of  the  machine  being  used,  the 


9 


forced  compatability  conditions  appear  to  be  non-zero.  The 


forced  conditions  at  O  ■  90°  and  270°  are  significantly  less 


than  the  rest  of  the  solution  and  therefore  satisfactory. 
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CHAPTER  6 
CONCLUSIONS 

The  apparent  simplicity  of  the  toroidal  shell  belles  the 
complexity  of  the  mathematics  to  analyse  this  shell 
structure.  The  analysis  of  simpler  shells  Cplates, 
cylinders,  and  spheres)  has  the  advantage  of  the  geometric 
curvature  remaining  either  constant  or  at  lease  maintaining 
its  sign  positive  or  negative.  The  geometric  curvature  of 
the  toroid  changes  from  positive  to  negative  as  ©  goes  from  0 
to  2n. 

The  impact  in  the  change  of  sign  in  the  radius 
of  curvature  about  the  Z  axis  has  been  seen  in  the  results 
of  this  work  and  in  applying  other's  solutions.  This  impact 
is  emperical  and  could  be  the  subject  for  further 
investigation.  The  first  blatant  statement  was  found  in  Tsui 
[241  where  he  states,  "...the  relevent  differential  equations 
to  assume  singular  solutions,  and  consequently,  problems 
involving  the  crowns  such  as  a  complete  toroid  cannot  be 


solved. "  The  culpret  was  writing  the  differential  equation 


in  terms  of  c Tsui's  nomenclature): 

r  *  b  <1  +  — — > 

2  sine 

(This  papers  nomenclature): 

p  =  r  (1  +  C^i/sin©) 

To  avoid  the  — - term,  Tsui  shifts  variables  -  both 

sine 

dependant  and  independant  -  and  provides  influence 
coefficient  matrices.  (These  are  outputs  from  numerical 
analysis. ) 

As  mentioned  earlier,  Tsui  addresses  the  significant 
differences  between  the  partial  toroidal  shape  and  the 
complete  toroid.  Others  may  have  addressed  partial  toroidal 
shapes  tangentially  but  Tsui  once  again  is  very  blunt  in  his 
statements.  Again  remember  Tsui  is  intending  his  work  as  a 
practical  handbook.  His  intended  audiance  wants  practical 
answers  to  real  problems.  Others,  Flugge  [25] 
specifically,  attempt  to  give  solutions  but  he  adds,  "...this 
solution  cannot  be  realized. .. because  it  again  leads  to  an 
incompatability  of  deformations. . . 

CURVATURE 

In  attempting  to  point  the  finger  at  why  the  solution 
is  so  difficult,  the  answer  comes  back  to  curvature.  This 
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work  has  seperated  curvature  into  geometric  (.unloaded)  and 

loaded  curvature.  This  writer  feels  that  the  geometric 

curvature  is  the  culpret.  Notice  the  expression  for  the 

radius  of  curvature  in  the  <f>  direction  given  in  chapter  one: 

p  -  — —  +  r 
sin6 

For  the  range  0°  <  ©  <  180°;  sin©  >  0  and  p  >  0.  Conversely, 
for  180°  <  0  <  360°;  sine  <  0  and  p  <  0.  At  ©  «  0°  and  180° 
sine  =  0  and  p  **  ±  Two  things  were  observed. 

A.  In  previous  solutions  obtained  in  this  study-  ie. 

worse  than  the  one  presented  -  the  stress  outputs  for 
the  two  regions  differed  markedly  between  0°  <  O  < 
180°  and  180°  <  ©  <  360°.  This  could  have 

been  due  to  a  lot,  of  reasons,  but  this  was  also 
coincident  with  . . . 

B.  Flugge’s  solution  for  the  u  displacement  t26J 
was  being  investigated.  This  tended  to  blow  up  from 
©  *  180°  to  360°/0°  due  to  a  term  ln<tan^).  Other 
problem  terms  existed  <cot©>  which  agrivated  the 
solution  at  the  points  6  ■  0°  and  180°.  In  just 
"plaving"  with  the  solution,  the  points  of  curvature 
change  kept  comming  back  as  critical  points  requiring 


further  attention. 

Something  happened  at  6  «  0°  and  180°.  In  a  cylinder 
these  two  points  are  no  trouble  -  but  the  radii  of  curvature 
are  «.  and  r  -  and  are  constant!  In  a  sphere  these  two  points 
are  no  problem  -  but  the  radius  of  curvature  is  r  in  all 
directions  -  and  constant!  In  the  toroid  the  radii  of 
curvature  at  both  points  are  p  *  ±  oo  and  r.  The  key 
difference  is  that  p  is  changing  from  p  >  0  to  p  <  0  and  at 
©  *=  0°  and  180°,  x^  (or  just  happens  to  be  passing  through 
zero. 

It  was  mentioned  earlier  in  Chapter  4  that  it  required 
two  reference  points  to  define  the  toroid.  Some  individuals 
may  take  exception  to  that  statement,  but  here  is  where  it 
comes  into  play.  Another  way  to  look  at  p,  the  radius  of 
curvature  in  the  <j!>  direction,  is  that  p  is  the  length  of  the 
arm  joining  the  point  in  question  (say  point  P>  and  the  axis 
of  symmetry  (the  Z  axis  in  our  case)  and  coincident  with  the 
local  unit  normal  vector  at  point  P.  The  angle  ©  can  be 
measured  either  at  the  axis  of  rotation  of  the  small  circle 
forming  the  shell,  or  at  the  intersection  of  the  line  p  and 
the  axis  of  symmetry  for  the  entire  toroid.  Now  when  ©  ■  50° 


lor  instance,  the  geometry  looks  like  that  shown  "in  Figure 
16.  As  6  decreases  toward  0,  tRe  end  of  p  must  travel 
farther  and  farther  down  the  axis  of  symmetry.  All  this  time 
p  is  getting  larger  and  x^  is  getting  smaller.  Vhen  ©  gets 
to  zero  plus  CO*)  the  situation  is  one  of  approximations.  p 
is  the  hypotenuse  of  the  right  triangle  with  the  right-angle 

p 

at  the  origin.  Sin©  *  -  from  simple  trigonometry,  but  as  p 
goes  to  sin©  goes  to  zero  and  hence  ©  must  go  to  zero 
also. 

Continuing  ©  around  to  0",  the  closest  point  on  the  axis 
of  symmetry  in  line  with  the  normal  vector  is  now  at  +«. 

This  now  gives  p  a  direction  opposite  to  that  of  the  local 
normal  vector,  hence  a  minus  sign,  and  as  ©  goes  further 
negative,  p  is  now  smaller  in  magnitude.  CSimilar  arguments 
can  be  made  for  the  bottom  crown  where  ©  *  180°>  The  sole 
reason  for  the  ±  <x>  trip  in  p,  and  hence  the  sign  change  in 
x^,  is  the  offset  R.  Somehow  this  must  be  related  to  the 
difficulty  in  solving  the  complete  toroid. 

GEOMETRY 

The  other  subtility  of  the  toroid  is  the  way  the 
mathematical  representat ion  of  the  surface  impacted  on  the 
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attempt,  to  integrate  for  energy  over  the  toroid  surface. 

This  work  specifically  avoided  some  problems  encountered  by 

other  investigators,  significantly  the  expression  for  the 

curvature  in  the  $  direction.  The  early  development  in  this 

work  specifically  avoided  the  term  C — - — ),  selecting  instead 

sin6 

an  expression  containing  the  term  Ca  +  sin©).  Force  balance 
and  constitutive  equations  were  developed  using  the  latter 
expression.  When  the  constitutive  equations  were  plugged 
into  the  energy  expression  <V>,  no  singularities  existed. 

This  was  good  news. 

Unfortunately,  the  Ca  +  sin©)  term  also  appeared  in  the 
expression  for  the  di f ferential  surface  area  CdS).  The  bad 
news  was  that,  many  of  the  elements  in  V  contained 
combinations  of  trigonometric  functions  not  readily  found  in 
the  CRC  Math  Tables  (271.  This  is  what  forced  the  numerical 
integration  ol  the  elements  in  V. 

This  heavy  dependence  of  the  trigonometric  functions  in 
describing  the  problem  has  an  astounding  impact  on  the 
selection  of  an  assumed  functional  for  the  energy  method 
solution.  Anv  combination  of  CA  ©  sin©)  or  C©A  sin©)  goes  to 
zero  when  integrating  from  0  to  2 n.  Remember  that  the  dS 
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term  contains  a  sin©  term  so  all  integrals  contain  at  least 
one  sin©  term. 

It  should  be  pointed  out  that  the  solution  to  a 
symmetric  problem  can  be  obtained  by  solving  only  one  of  the 
unique  parts  of  the  problem.  This  was  not  done  in  this  case 
because  it  was  hoped  that  the  program  being  developed  would 
be  more  general  in  nature,  actually  leading  to  orthotropic 
analysis  of  a  stiffened  toroid. 

The  only  redeeming  quality  of  the  integrals  also  comes 
from  geometry.  Geometry  terms  containing  cos©,  cos2©,  sin©, 
sin2©,  and  <a  +  sin©;  and  thler  products  appear  in  numerators 
and  denominators  through  out  the  integrals.  Mister  Simpson 
and  the  digital  computer  were  allowed  to  seperate  the  wheat 
from  the  chaff. 

An  attempt  was  made  to  develope  a  non-trigonometric 
functional  which  would  satisfy  the  conditions  of 

connectivity  stated  earlier  and  to  meet  connectivity  at 
©  *  0°/360°.  The  conditions  to  be  met  were: 


f<0>  -  f (360) 
f’cO)  ■  f ' (360) 


f (90)  «  0 
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r  C90)  ■  o 


f(270>  -  0 
f’<270>  -  0 

This  resulted  in  eight  conditions.  A  functional  of  the 
following  form  was  generated: 

f <  m0)  «  A  +  EKm0>  +  C(m0>2  +  DCm©)3 

+  E(m©)4  +  FCm©)5  +  G(m©>6  +  HCm©)7 
A  was  arbitrarily  set  to  1  (to  force  a  non  zero  solution  at 
0  «  0°).  Unfortunately,  by  the  time  f(mq)  was  processed, 
the  functional  looked  the  same  for  all  values  of  m.  The 
coefficients  all  changed  as  m  changed,  but  the  shape  of  the 
fuhctional  remained  unchanged  from  one  value  of  m  to  the 
next  because  of  the  eight  conditions  imposed  upon  the 
functional.  This  lost  the  variety  of  the  solution  that  one 


function  would  satisfy  the  form  of  the  differential  equation 
for  the  toroid  and  hence  yield  an  exact  solution.  <It  was 
partially  the  anticipation  of  using  such  a  complex  function 
that  drove  the  layout  of  the  program  (see  Appendix!)  to  have  a 
seperate  generating  section  for  the  functional.?  The 
Legendre  function,  or  any  other  hypergeometric  function,  was 
not  used  due  to  compatability  conditions  being  harder  to 
achieve.  The  next  step  in  this  progression  would  be  the 
assumption  of  a  combined  trigonometric  function  for  u  and  w, 
Cie.  A  sinCmO)  +  B  cos<m©)>. 
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Figure  17 

Generated  Functional 
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CLOSING 


In  Chapter  1  it  was  stated  that  the  geometric  curvature 

and  the  impact  of  shell  geometry  on  the  assumed  functional 
were  the  two  problem  areas  in  this  thesis.  Of  the  two,  the 
former  is  the  more  significant  in  this  writers  opinion.  The 
analysis  behind  Figure  15  is  crude,  but  it  is  an  attempt  to 
represent  in  graphics  and  mathematics  what  the  shell  is 
doing.  Being  unable  to  do  this  is  what  makes  the  solution  so 
difficult. 

In  a  dynamics  class,  one  student  asked  the  professor  if 
another  mathematical  technique  could  be  used  to  analyse  a 
spring.  The  professor  thought  briefly,  then  answered,  "You 
can  do  all  the  math  you  want  but  you  have  to  think  like  a 
spring.  The  spring  knows  what  it  is  doing. "  This  writer 
feels  a  little  short  on  thinking  like  a  toroid,  but  also 


feels  in  good  company. 


APPENDIX 


THE  PROGRAM 

This  program  is  to  set  up  the  geometry,  the  energy 
expression.  Simpson’s  multipliers,  and  perform  the  matrix 
manipulations  to  solve  for  the  displacements  <u>  and  <w>. 


50  CLS: PRINT: PRINT: INPUT  "WHAT  OUTPUT  FILE  NAME  ";A$ 

60  CLS 

1000  NN=  50  *This  is  the  number  of  steps  to  be  used  in  the 
Simpson's  integral 

1010  MM=1 0  *This  is  the  number  of  terms  in  the  series  for 
displacements  u  ans  w. 

1020  NU=.3  *Poisson's  ratio 

1030  R=2  *R  from  the  text. 

1040  A=4  *a  from  the  text. 

1050  H=.02  *t  from  the  text. 

1060  E=3E+0?  #Young’s  modulus 
1070  P=-100  »Kpressure  in  psi. 

1080  CLS: PRINT: PRINT: PRINT 

1090  PRINT  "The  following  are  the  default  values  in  the 


program: 

1100  PRINT: PRINT 
1110  PRINT  "GEOMETRY" 

1120  PRINT  "Radius  of  rotation 

1130  PRINT  "Radius  of  the  circle 

1140  PRINT  "Thickness  of  the  shell 

1150  PRINT: PRINT 

1160  PRINT  "MATERIAL- 

1170  PRINT  "Poisson  s  ratio 

1180  PRINT  "Modulous  of  elasticity 

1190  PRINT: PRINT 

1200  PRINT  "SOLUTION" 

1210  PRINT  "Number  of  steps  in  the  Simpsons" 
1220  PRINT  "  integral 

1230  PRINT  "Number  of  terms  to  be  used  in" 
1240  PRINT  "  the  series  approximation 
1250  PRINT  "Pressure 


"  inches" 

inches" 

inches" 


psi 


<  i 


1260  PRINT: PRINT  „ 

1270  INPUT  "Do  you  wish  to  change  any  of  these  CY/N):  :S8 
1280  IF  SS="N"  THEN  1520 


1290  CLS 

1300  PRINT: PRINT 
1310  PRINT  "GEOMETRY" 

1320  INPUT  "Radius  of  rotation  Cinches)  :";RR 

1330  INPUT  "Radius  of  the  circle  Cinches)  :";R 

1340  A=RR/R 

1350  IF  A  >  1  THEN  1380 

1360  CLS: PRINT"RADIUS  OF  ROTATION  MUST  BE  GREATER  THAN  THE 
RADIUS" 

1370  PRINT  "OF  THE  CIRCLE. ": GOTO  1300 

1380  INPUT  "Thickness  of  the  shell  Cinches)  :";H 

1390  PRINT: PRINT 


1400  PRINT  "MATERIAL- 

1410  INPUT  "Poisson  s  ratio  :";NU 

1420  INPUT  "Modulus  of  elasticity  :":E 

1430  PRINT: PRINT1440  PRINT  "SOLUTION" 

1450  PRINT  "Number  of  steps  in  the  Simpsons" 
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1460  INPUT  "  integral  :";NN 

1470  PRIM  "Number  of  terms  to  be  used  in" 

1480  INPUT  "  the  series  approximation  :";MM 
1490  INPUT  "Pressure  (psi)  :‘‘;P 

1500  ' 

1510  ' 

1520  'THIS  SECTION  IS  TO  SET  UP  THE  PARAMETERS  FOR  THE 
INTEGRATION  SECTION 

1530  ’ 

1540  DT=6. 283185307^/NN 
1550  1=NN+1 

1560  DIM  PC  1.3, MM)  ^Defines  the  functions  <p. 

1570  DIM  JCI)  ♦Simpson’s  multiplier 

Arrays  MD.  ME.  MF,  and  MG  correspond  to  the  matrixes  in 
the  partial  of  V  with  respect  to  uk  and  *k.  MP  is  the  matrix 

representing  the  energy  due  to  the  pressure  term.  MS  is  a 
scratch  matrix.  M*~  is  a  series  of  matrixes  used  in  the 
matrix  manipulation  subroutines.  MV  and  MU  are  the  answers. 

1580  DIM  MDCMM.MM): DIM  MECMM. MM) : DIM  MFCMM, MM) 

1590  DIM  MGCMM.MM): DIM  MSCMM, MM) 

1600  DIM  MX^CMM. MM) : DIM  MY*CMM, MM) : DIM  MZ^CMM, MM) 

: DIM  MI^CMM. MM) 

1610  DIM  MPCMM) : DIM  MVCMM) : DIM  MUCMM) 

1620  CLS: PRINT: PRINT: PRINT  "GENERATING  FUNCTIONS" 

1630  FOR  I  =  1  TO  NN+1 
1640  T=C 1-1 )*DT 
1650  FOR  M  =  1  TO  MM 

1660  PCI.1 .M)=SINCM*T)  ’This  is  the  function 

1670  PCI.2,M)=M*C0SCM*T)  'This  is  the  first 

derivative 

1680  PCI  ,3.N)=-M'2*SINCM*T)  'This  is  the  second 

derivative 

1690  NEXT  M 
1700  KI)=4 

1710  IF  C 1/2-INTC 1/2) ) >0  THEN  JC I )=2 

1720  NEXT  I 

1730  JC1 )=1 : JCNN+i )=1 

Integrate  using  Simpson’s  Rule 

2000  ’INTEGRATE  FROM  0  TO  2*PI 
2010  IC=3 .  14159*CE*H/Cl-NU'2))*DT/3 
2020  BC=H'2/C6*R'3) 

2030  PC=-2*3 .141 59*P*R'  2*DT/3 
2040  FOR  I  =  1  TO  NN+1 
2050  T=C 1-1 )*DT 
2060  S=SINCT) 

2070  C=C0SCT) 

2080  CLS: PRINT: PRINT: PRINT-INTEGRATING  ON  STEP  ":I 
2090  PRINT  "  ANGLE  " : T*57 . 296 ; "  DEGREES" 

2100  AA=A+S2102  HH=H' 2/C12*R^2) 

2104  CC=C/AA 
2106  CA=2*NU*CC 
2108  CB=CC' 2 
2110  CD=2*NU*S/AA 
2112  AX=2*S*CVAA'  2 
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K 


B 

I 

<•: 

JS 

Ei 

r 


I 


2130 

2140 

2150 

2160 

2170 

2180 

2190 


2200 


2210 


2220 


i 

\\ 

2230 

NEXT  M 

V~ 

2240 

NEXT  k 

l  < 

2250 

NEXT  I 

& 

2260 

* 

2270 

y 

i 

2280 

FOR  k  =  1  TO  MM 

P 

2290 

MF(k)-FC*MP(k) 

2300 

FOR  M  =  1  TO  MM 

2310 

MDCk, M)=IC*MDCk,  M) 

2320 

MECk, M)=IC*MECk, M) 

2330 

MFC  k .  M)  =  IC*MFC  k .  M) 

2340 

MGCk.M>=IOMGCk.M.> 

4.-  ■ 

2350 

NEXT  M 

T' 

2360 

NEXT  k 

3000 

'SOLVE  THE  MATRICES 

FOR  K  =  1  TO  MM 
MPCk)=MPCk)-JCI)*PCI.l,k)*AA 
SQ=1/AA'  2 
TR=2*NU/AA 

BX=H'  3*0  +2*NU*C/AA+C  C/AA) '  2  )  /C 1 2*R'  2 ) 

FOR  M  =  1  TO  MM 

MDCk.M)=MDCk,M)+JCI)*CCl+HH)*(2*PCI.3,R>*PCI,3,M) 
+CA*PC  1.3,  k)*PC  I , 2 , M$ 

+CB*2*P( I ,2, k)*P(I, 2 ,M)))*AA 
MECk,  M)=MECR,  M)+JC I)*(2*Pc 1 . 1 ,M)*P(1 , 3,  k) 
+CA*CPCI,  3,k5*Pa.l,M5*S 
+P(I,l,M3*f>(I,2,k}*fc:> 
+AX*f>cL2.k)*PCi,l.M)  _ 

-HH*C2*Pct ,3, M)*P(1 , 3, k) 

+CA*CPd.3.k)*P(5:.2,M) 

+PCI,3,M5*P(I,2.k>) 

+2*C§*Pc  1 , 2 .  k5*PC  1 . 2 .  M)  )  >*AA 
MFC k .  M )  =MFC k .  M )  +  J C I )*( 2*PC 1 , 3 ,  M) *PC  I  ,  1 .  k> 
+2*NU*CP(  1 . 1 .  k;*P(i ,3 .  M$*s 
+PC I , 2 . M)*PC I ,1 ,k)*C)/AA 
+AX*Pc:r  2,M)*P(I,l,k) 
-HH*C2*Pct.3,M)*P(i ,3.k) 

+CA*(P(i . 3,M)*P(I . 2,  k) 

+PC I . 3,k5*Pci  2, Mi) 

+2  *CB*P  Cl , 2 . M  5  *f>C 1 , 2 . k> ) ) *AA 
MGCk.M)=MGCk,M)+J(i:>*<2*Pa.l,M>*Pa.l,k) 
+2*CD*P( I, 1 ,  ff)*PC  I  1,D 
+2*PC I , 1 . fc)*PC I . 1 , fo*S~2/AA~ 2 
+HH*C2*PC  I  3, M)*PC  1 . 3,  k> 

+CA*Cpct.3,M)*Pci.2,k> 
+PCI.3.k$*Pci.2.M$) 

+2*CB*PC  1 , 2 ,  M5*P(  1 , 2 ,  k) ) )*AA 


3010  CLS: PRINT: PRINT: PRINT  "SOLVING" 

3020  ’ 

3030  FOR  k  =  1  TO  MM: FOR  M  =  1  TO  MM: MX#(k, M)=MDCk, M) : NEXT 
M: NEXT  k 

3040  G0SUB  10000  ’  *****  Invert  MD  **** 

3050  FOR  k  =  1  TO  MM: FOR  M  =  1  TO  MM 
3060  MY*Ck. M)=MECk. M) : MX*Ck, M)=-MI^(k, M) 

3070  NEXT  M: NEXT  k 

3080  G0SUB  11000  '  ****  C-MD'-1)*ME  **** 

3090  FOR  k  =  1  TO  MM: FOR  M  =  1  TO  MM 
3100  MY*Ck.M)=MZ*Ck.M>: MX*Ck. M)=MFCk, M) 

3110  MS'Ck.M)=MZ^Ck.M)  ’  ****  Saving  this  for  Um 

solution  **** 

3120  NEXT  M : NEXT  k 
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3130  GOSUB  11000  '  ****  MF*(-MD' -1 )*ME  **** 

3140  FOR  k  =  1  TO  MM: FOR  M  =  1  TO  MM 
3150  MIXk,  M)=MG(k,  M) :  MX*(k, M>=MZ#Ck, M> 

3160  NEXT  M: NEXT  k 

3170  GOSUB  12000  ’  ****  MG  -  MF*(-MD~-1 >*ME  **** 

3180  FOR  k  =  1  TO  MM: FOR  M  =  1  TO  MM: MX*(k. M)=MZ*(k, M) 

.-NEXT  M:  NEXT  k 

3190  ’ 

3200  ’ROV  REDUCTION 

3210  ’Input  is  MX^Ck.M)  and  MPCK> 

3220  ’Output  is  the  altered  MX^CK.M)  and  MPCK) 

3230  FOR  S  =  1  TO  MM-1 
3240  FOR  k  =  S+l  TO  MM 
3250  Q=-MX*(k.S)/MX*CS.S) 

3260  MP(k)=Q*MP(S)+MP(K) 

3270  FOR  M  =  1  TO  MM 

3280  MX^(k.M)=0*MX^CS,M)+MX#Ck.M) 

3290  NEXT  M 
3300  NEXT  k 
3310  NEXT  S 
3320  ’ 

3330  ’SOLVE  FOR  MVCM) 

3340  ’Continuing  on  using  MX^(k.M)  and  MP(k)  from  before 
3350  FOR  S  =  MM  TO  1  STEP  -1 
3360  MV(S)=MP( S) 

3370  FOR  M  =  S+l  TO  MM 

3380  MVT(S)=MV(S)-MV<M)*MX*(S.  M) 

3390  NEXT  M 

3400  MV/(S)=MV(  S)/NX*t S .  S) 

3410  NEXT  S 
3420  ; 

3430  ’  SOLVE  FOR  MUCk) 

3440  FOR  k  =  1  TO  MM 

3450  FOR  M  =  1  TO  MM 

3460  MUC  k  )=MU(.  k) +MSC  k  .  M  ) *M  V(  M  > 

3470  NEXT  M 
3480  NEXT  k 

3490  ’  PRINT  OUT  RESULTS 
3500  PRINT  Vm  Um  " 

3510  FOR  M  *  1' TO  MM 

3520  PRINT  M . MVCM) . MHf M> : PRINT 

3530  NEXT  M 

3540  OPEN  “O" . #1 . AS  ‘•‘This  section  puts  the  output  onto 

disk . 

3550  PRINT  P\  ,  "NN  =  " ;  NN  :  "  MM  =  MM 
3555  PRINT  *1 . "RADIUS  =":A*R;"  radius  =";R 

3560  PRINT  Pi/'  ” 

3570  PRINT  *1  ."*"/’  Vm  " . "  Um  " 

3580  FOR  M  =  1  TO  MM 

3590  PRINT  PI .  M.MV(M),MU(M): PRINT 

3600  NEXT  N 

3610  CLOSE 

3620  END 
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This  section  contains  routines  used  in  the  main  program. 

10000  ’ INVERSION  OF  A  MATRIX 

10010  'Input  is 

10020  'Output  is  MI*Ck,M) 

10030  FOR  k  =  1  TO  MM 

10040  MI#(k.k)=l 

10050  NEXT  k 

10060  FOR  S  =  1  TO  MM-1 

10070  FOR  k  =  S+l  TO  MM 

10080  Q=-MX#(k.S)/MX*(S.S) 

10090  FOR  M  =  1  TO  MM 

10100  MX#(k,  M)=Q*MX>?(S.  M)+MX>*(k,  M) 

10110  MI^Ck.M)=Q*MI^CS.M>+MI#(k.M> 

10120  NEXT  M 
10130  NEXT  k 
10140  NEXT  S 

10150  FOR  S  =  MM  TO  1  STEP  -1 
10160  FOR  k  =  S-i  TO  1  STEP  -1 
10170  Q=-MX*(k.S>/MX#CS,S) 

10180  FOR  M=MM  TO  1  STEP  -1 
10190  MX*(k.M)=Q*MX*(S,M:>+MX#(k.M) 

10200  Ml^Ck. M)=0*MI^CS. M)+MI#(k. M) 

10210  NEXT  N 

10220  NEXT  k 

10230  NEXT  S 

10240  FOR  k  =  1  TO  MM 

10250  FOR  M  =  1  TO  MM 

10260  MI^(k.M)=MI^(k.M)/MX^Ck.k> 

10270  NEXT  M 
10280  NEXT  k 
10290  RETURN 

11000  MULTIPLY  TVO  MATRICES 

11010  'Multiply  MX^Ck.M)  times  Ml>V(k,M)  in  that  order! 

11020  'The  output  is  MZ^lk.Mi 
11030  FOR  k  =  1  TO  MM 

11040  FOR  M  =  1  TO  MM11050  MZ*(k.M>=0 

11060  FOR  S  =  1  TO  MM 

11070  MZ^Ck.  M)=MZ*(k,  M)+MX«‘(k.  S^MY^CS,  M) 

11075  IF  ABS(MZ^(k,M)X.  0000005  THEN  M^(k.M)=0 

11080  NEXT  S 

11090  NEXT  M 

11100  NEXT  k 

11110  RETURN 

12000  'ADD  TVO  MATRICIES 
12010  'Add  MXtfCk.M)  to  MY*(k.M) 

12020  'Output  is  MZ*(k,M) 

12030  FOR  k  =  1  TO  MM 

12040  FOR  M  =  1  TO  MM 

12050  MZ*Ck.M)=MX*(k,M)+MY*(k,M) 

12060  NEXT  M 
12070  NEXT  k 
12080  RETURN 
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